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Abstract 

We introduce an approximation scheme for the Hodgkin-Huxley 
model of nerve conductance which allows to calculate both the speed 
of the traveling pulses and their shape in quantitative agreement with 
the solutions of the model. We demonstrate that the reduced problem 
for the front of the traveling pulse admits a unique solution. We obtain 
an explicit analytical expression for the speed of the pulses which is 
valid with good accuracy in a wide range of the parameters. 
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1 Introduction 

Understanding the mechanisms of the propagation of the nerve activity is one 
of the fundamental problems in biophysics. The simplest example of such a 
propagation is a single solitary traveling pulse of action potential in an axon 
(|Katz, 1966|) . Today it is well established that the changes of the membrane 
potential in nerve tissue are the result of the complex dynamics of the ionic 
currents through voltage-sensitive channels ( [Katz, 1966| ). The first detailed 
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quantitative measurements of the ionic currents were performed by Hodgkin 
and Huxley in the early 50-s (Podgkin and Huxley, 19521) . Using the voltage 
clamp technique, they were able to measure the kinetics of Na"*" and K"*" cur- 
rents in the squid giant axon. This led them to a set of differential equations 
which describe the dynamics of the action potential. Furthermore, by com- 
bining these equations with the cable equation for spreading of current in the 
axon they were able to calculate the shape and velocity of the propagating 
action potentials ( Podgkin and Huxley, 195^ ; [Huxley, 1959|) . The predictions 
of their model turned out to be in a remarkably good agreement with the 
experimental observations. 

The reason that the model introduced by Hodgkin and Huxley (the HH 
model) admits quantitative comparisons with the experiments is that it con- 
tains detailed information about the volt age- dependent kinetics of the Na"*" 
and K"*" channels. Naturally, this makes the models quite complex and in- 
tractable analytically. So far, the basic tool for studying the HH model were 
numerical simulations. Note that because of its complexity, it was not until 
recently, with the advent of very fast computers, that the simulations of the 
HH model could be done routinely. Even then, one is still required to do sim- 
ulations for each set of the parameters of the model. Therefore, analytical 
studies which would give functional dependences of the main parameters of 
the action potentials on the parameters of the model are still highly desirable. 

The early analytical works on the HH model relied on the strong sep- 
aration of the time scales of the (fast) activation and (slow) inactivation 
processes. These studies made an assumption that the Na^ activation is the 
fastest process and can be eliminated adiabatically, what amounts to assum- 
ing that the sodium activation variable m = rrioolV), where moo{V) is the 
resting value of m at a given membrane voltage V ( [FitzHugh, 1961| ; [Casten 



et al., 1975| ; [Carpenter, 1977| ; [Carpenter, 1979[ ). This leads to a cubic-like 
nonlinearity in the equation for the membrane potential. By further assum- 
ing that the Na"*" inactivation and K"*" dynamics are much slower than the 
Na"*" activation, the problem of the action potential propagation reduces to a 
single reaction-diffusion equation for the front of the action potential ( |Uas-[ 
ten et al., 1975 ). A number of simpler models (FitzHugh-Nagumo type) 
with similar properties had been introduced to mimic the behavior of the 
membrane ([FitzHugh, 1961[ ; [Nagumo et al., 1964[ ; [Rinzel and Keller, 1973[ ; 
[Casten et al., 1975| ; [Jones et al., 199l[ ). The latter in fact became quite 
popular for explaining traveling waves phenomena in a variety of excitable 
systems in physics, chemistry, and biology (jVasiliev et al., 1987[ ; [Murray 
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1989; llVlikhailov, [Kerner and Usipov, 199^ . 

Although this kind of analysis leads to a qualitative explanation of the 
excitability of the nerve membrane, it fails to give any quantitative predic- 
tions for the speed of the propagating action potentials. It also predicts that 
the traveling wave should have the form of a broad excitation region with 
the sharp front and back. This is in contrast to the observations in which 
the pulse is a narrow localized region of excitation (a spike). The reason for 
this is that in reality it is typically the membrane potential rather than the 
Na^ activation that is the fastest process. For example, in the squid giant 
axon the time constant of the membrane potential change is Ty ~ 0.01 ms, 
whereas the time constant of the Na^ activation is roughly ~ 0.2 ms. So, 
the FitzHugh-Nagumo type models are in fact not adequate for any quan- 
titative predictions of the characteristics of the action potential. Also, they 
only qualitatively reveal the mechanism of the wave propagation. 

In the present paper, we introduce an approximation scheme that does 
take into account this relationship between the time scales. We will construct 
an approximate solution for a single traveling pulse in the HH model that is 
in the quantitative agreement with the solutions of the full HH model. We 
will investigate the structure of the front of the traveling pulse and show 
that it is substantially different from the conventional case of the FitzHugh- 
Nagumo type models. We will also obtain an explicit analytical expression 
for the speed of the pulses which agrees with the results of the simulations 
of the HH model within 20% accuracy in a wide parameter range. Using the 
obtained solutions, we will construct an approximate solution for the entire 
pulse that is also in quantitative agreement with the solutions of the full HH 
model. 



2 The Hodgkin-Huxley model 

In the following, we will use the version of the HH model which was originally 
introduced by Hodgkin and Huxley to study the behavior of the squid giant 
axon ( |Hodgkin and Huxley, 19521) and later adopted by many researchers as 
a benchmark of the models of nerve activity. Namely, we will consider the 
following equations 

dV a d'^V 

C— = —j-^ + g^.rrv'h{y^^-V)+gKn\VK-V)+giiyi-V),(l) 
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amiV)il-m)-/3miV)m, (2) 
ah{V){l - h) - f3h{V)h, (3) 
an{V){l-n)-(3n{V)n. (4) 



Here, as usual, Eq. (P is the cable equation for the membrane potential V, 
with C = 1 /iF/cm^ the membrane capacitance per unit = 238 /im is 

the radius of the axon and p = 35.4 f2 x cm the resistivity of the intracellular 
space; (^Na = 120 mf2~^/cm^, = 36 mf2~^/cm^ are the conductances of the 
open Na"*" and K"*" channels per unit area; Vne = 115 mV and Vk = —12 mV 
are the equilibrium potentials of Na"*" and K"^, and gi = 0.3 mfi^^/cm^ and 
Vi = 10.5989 mV are the leakage conductance per unit area and the leakage 
voltage, respectively. With these definitions the resting potential K = 0. 
Similarly, m and h are the activation and the inactivation variables for the 
Na"*" channels, respectively; n is the K"*" inactivation variable; the rates am,h,n 
and Prn,h,n as functions of V at temperature T = 6.3° C can be found in 
( Hodgkin and Huxley, 1952|) (note that ( |Hodgkin and Huxley, 1952 ) have an 



opposite sign convention for V). The temperature changes are accounted for 
by a factor = 3(r-6-3)/io multiplying all a's and /3's, T is in degrees Celsius. 
The lengths are measured in centimeters and the times in milliseconds. The 
voltage V is measured in millivolts. 

Equations (0) - (H) constitute a closed system of partial differential equa- 
tions that quantitatively describe the changes in the membrane as functions 
of time and space. Let us emphasize that their ingredients are obtained by 
measurements and fitting of the parameters to the actual experiments. So, it 
is important to understand the relationships between the characteristic pa- 
rameters, namely the time scales, in this system. From the functional form of 
a's and /?'s ( [Hodgkin and Huxley, 1952D we can make the following estimates 
for the time constants Tm,h,n for m, h, and n, respectively, at T = 6.3° C: 

Tm ~ 0.2 msec, (5) 
~ 5 msec, (6) 
r„ ~ 3 msec. (7) 

Also, from Eq. (|l|) one gets the following estimate for the time scale ry of 
the variation of the voltage, assuming that all the Na"*" channels are open: 

Ty ~ C/(7Na ~ 0.01 mSCC. (8) 
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One can see that the following hierarchy of time scales holds in the system: 



Tv < Tv„ < Th, r„. (9) 

The first inequality holds better for sufficiently low temperatures, and re- 
mains qualitatively correct up to the temperatures T ~ 30° C at which the 



pulses fail to propagate in the HH model (puxley, 19591) . As we already 
pointed out in the Introduction, this is an important property of the system 
which is not taken into account in most of the analytical studies of the HH 
model. In the following, we will use this hierarchy of time scales to introduce 
the approximation scheme for the traveling pulses in this model. 

Another important point about the HH model is the fact that the very 
nonlinearities in Eq. ([^), namely, the powers with which the variables m, h, 
and n enter the equation are determined experimentally ( [Hodgkin and Hu5c^ 
ley, 1952|) . Furthermore, these powers correspond to the number of particles 
involved in the operation of the respective channels and therefore represent 
significant physical quantities. As will be seen below, these powers in fact 
play crucial role in our studies. 

Before going to the analysis of the traveling pulses, let us discuss the 
basic physics of the excitability in the HH model. In the rest state, the Na"*" 
channels are basically closed, at \^ = the equilibrium values for m and h are 
mo — 0.05 and ~ 0.60, respectively, while the channels are partially 
open, with no — 0.32. If, by applying an external stimulus, the membrane 
voltage V is increased to ~ 10 mV, the Na"*" channels will start opening on 
the time scale of order r^. The infiux of the Na^ ions will in turn lead to 
the increase of the membrane potential V on the time scale intermediate 
between Ty and (see below), resulting in the positive feedback. The 
membrane potential V will come close to the resting potential V^^, while 
the Na"*" channels will become mostly open with m ~ 1. During this time, 
the slow inactivation variables h and n will remain almost unchanged. After 
that, the slow inactivation variables h and n will start to react, closing the 
Na^ and opening the K+ channels, which will drive the potential V back 
to equilibrium. In the spatially extended system the diffusive spreading of 
the current in front of the excitation region in the axon will provide the 
sustaining force for the propagation of the pulse along the axon. In that 
sense, from the physical point of view the traveling pulse in the nerve axon is 
a classical example of an autosoliton — self-sustained solitary inhomogeneous 
state in an active dissipative system whose existence is determined only by 
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the nonlinearities of the system and not the initial conditions ( [Kerner and 
Osipov, 1991) . 



3 Solitary pulse 

We are now going to construct an approximate travehng wave solution in the 
form of a solitary pulse, using the ideas introduced in the preceding section. 
Let us introduce a self-similar variable z = x — ct, where c is the propagation 
speed of the pulse. Then, Eqs. (|l|) - @) for a traveling wave with speed c 
will become 



2p dz^ dz 

(10) 

QTn. 

+ araiV){l - m) - f3m{V)m = 0, (11) 
+ aH{V){l -h)- h{y)h = 0, (12) 

djTi 

+ «„(V^)(1 -n)- /3niV)n = 0. (13) 

The boundary conditions for these equations are 

V{±oo) = 0, m(±oo) = mo, h{±oo) = Hq, n(±oo) = no, (14) 

where mo, h^, and uq are the values of m, n, and h in the rest state, respec- 
tively. For the chosen functions a and jS the rest state ^ = is unique and 
stable. 

The solution of Eqs. ( [TI1| ) - (|T^) in the form of a traveling solitary pulse 
obtained numerically at the "standard" temperature T = 6.3° C is shown in 
Fig. |I[ From this figure one can see several features of the solution which we 
will use in the approximation scheme that we are going to construct. First, 
observe that the length scale of the rise of the potential is substantially 
smaller than that of the fall of the potential. Second, during the rise of the 
potential the variables h and n remain almost unchanged at their resting 
values ho and Uq. Third, in front of the spike the value of m (that is, mo) is 
practically zero. 

Let us use the above facts to simplify Eqs. (|10|) - (|13|). Since the values 
of h and n change little in the front of the spike, we may replace them by 
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Figure 1: The numerical solution of Eqs. ([T0| ) - (13) at T = 6.3° C 



their values ho and tiq at rest and disregard Eqs (|T2D and (|^). Furthermore, 
since the value of mo is very small, with very good accuracy we may assume 
it to be zero. Therefore, in the rest state we will have 5'k^o^ + 5''^ — with 
very good accuracy, so these terms drop out of Eq. (|ToD . Also, the coefficient 
multiplying —V in the last two terms of Eq. (|TD|) is of order 0.7 and is much 
smaller than the contribution from the term g^arn?h during the rise of the 
potential when m is not close to zero, so these terms can be dropped as 
well. What we are then left with is a an equation for V coupled only to the 
equation for m with a number of terms dropped. Observe that since V is 
much faster than m when m ~ 1, the value of m has to be sufficiently small 
in order for the nontrivial collective dynamics of V and m to be possible. 
This allows to further simplify Eq. (0) by neglecting the terms proportional 
to m. After making all these approximations, we are left with the following 
set of equations 

« + + g^.m'hoiV^. -V)=Q, (15) 



2p dz^ dz 
dTfi 

c-^ + am{V) - a„(0) = 0, (16) 
dz 

instead of Eqs. (|l^) and (|ll|). Note that we added a term — am(0) to Eq. (|16D 
in order for this equation to be consistent with the approximate boundary 
conditions ahead of the spike front 

m(+oo) = 0, V{+oo) = 0, V;(+cx)) = 0, (17) 
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where = dV/dz. We can do this in our approximation scheme since the 
value of arn(O) is in practice very small compared to am(^a)- 

Let us assume that the characteristic value of m in the front is m <^ 1 
and the characteristic width of the front is I. Then, since all the terms in 
Eqs. ([TsD and ( |16|) should be of the same order of magnitude, we obtain the 



following estimates 

a cC , _o cm . . 

— — gNahom , — ~ am(\/Na)- (18) 

From these one can also estimate the characteristic time scale for the rise of 
the potential in the pulse as r ~ //c ~ (C'/a^fy'Na^o)^^'') where 

ttm = am(VNa) " am(0), (19) 

so 

3 \ 1/4 



r 



One can see from this equation that the dynamics in the front of the traveling 
pulse will indeed occur on the time scale intermediate between Ty and Tm- 

From the estimates above we immediately conclude that in the traveling 
spike 

a^ag^^h,\ m = ^) , I = -—J- -) ,(21) 



where c is a constant of order 1. Substituting the parameters of the HH 
model at T = 6.3° C, we see that rh ~ 0.6, what corresponds to the relevant 
quantity rh^ ~ 0.2, which is indeed rather small. 

Indeed, let us introduce the following new variables: 

z (?m V _, ,du 

4 = -, s = ^—, u=—, V = a{u)—, (22) 
I m l/Na ds 

where 

_/ N am(H^aM) - am(0) 

«(W) = 77^^ 77^- (23) 

am(VNa) - am(0) 

The latter is plotted in Fig. |^. Note that this way the time scale in Eq. (p!6[) 
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U 

Figure 2: The dependence a{u). 

has been absorbed into the constant a^. 

In terms of the variables introduced in Eq. (^) and after a few manipu- 
lations one can rewrite Eqs. (0) and ([iBD in the following form 

5('«)^ = !', (24) 



^ -i. (26) 

where now s is an independent variable. These transformations can be done 
for < M < 1 since a{u) is always positive for -u 7^ (see Fig. ^. Note that 
these equations do not have any ^ dependence in their right-hand side, so it 
suffices to solve Eqs. and (^) only. The solution for ^(s) can then be 
obtained by a simple integration. 

The problem now became substantially simpler because instead of solving 



the nonlinear boundary value problem for Eqs. (|T0|) - (|15), one now needs to 
solve the initial value problem for Eqs. (|2^) and (|25|). Indeed, according to 
Eq. (p!?!), when z +00 we have m — > 0, so s ^ as ^ — > +00. This means 
that M = and f = at s = 0, since also du/d^ = —cv — as ^ — >■ +00 
[see Eqs. (IT^), (P^), and (P^)]. One should though be careful to specify what 



exactly happens near s = since a(0) = 0. To do this, let us divide Eq. 
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by Eq. (mj. We get 



du (?v 

When s — > 0, we have dv/du 1, provided that v{s) has a non-zero deriva- 
tive at s = (the latter follows from the physical considerations). Therefore, 
according to Eqs. (p^) and (pS]), as s ^ 0, the solution will behave like 



where the prime means differentiation. Here we expanded a{u) in the neigh- 
borhood of zero and took into account that a'{0) ^ 0. 

It is not difficult to see from Eq. (^71) that for < m < 1 and v > the 
projection of the phase trajectory on the u — v plane will lie below the line 
u = V. Since du/ds > and there are no fixed points in this region of u 
and V, the solution u{s),v{s) will cross either the line m = 1 or f = 0. By 
direct inspection of Eqs. ( ^4]) and (|25| ) one can see that this intersection is 
transversal. Observe that the intersection of the lines u = 1 and w = is a 
fixed point in this plane. 

According to Eqs. (p^ and (^S]), once the solution leaves the region 
bounded by the lines u = v, u = 1 and f = 0, it can never come back 
to this region. Indeed, if the solution crosses the line m = 1 at some f > in 
the u — V plane, it will then have dv/ds > for all s, so v will stay positive. 
If on the other hand, the solution crosses the line f = at some u < 1, we 
will have both dv/ds < and du/ds < afterwards. In fact, it is easy to see 
that the only trajectory on which u and v will remain bounded for all s is the 
one which connects the point u = 0,v = with u = l,v = 0. Therefore, this 
trajectory is precisely the one that corresponds to the front of the traveling 
pulse. 

Of course, not all the speeds c will produce this kind of trajectory. It is 
clear that if c is very large, the trajectory will cross the line m = 1 in the u — v 
plane at v close to 1. On the other hand, if c is very small, the trajectory 
will cross the line f = at very small u. Figure |^ shows the results of the 
numerical solution of Eqs. (^) and ( P5[ ) at different values of c. From this 
numerical solution we found that the trajectory that connects u = 0,v = 
and u = l,v = 1 exists only for a unique value of c = c*. 

In fact, it is possible to prove that such a trajectory indeed exists and 
is unique at a unieque value of c. To do this, let us see what happens 
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Figure 3: The numerical solution of Eqs. ( ^4] ) and ( pS] ) in the u — v plane 
v{u) at different c. 

with the trajectory as the value of c is changed. For convenience, we will 
use u instead of s as an independent variable. Let vq{u) and sq{u) be a 
trajectory in the region bounded by m = f , m = 1 and f = with the initial 
conditions fo(0) = 0, so(0) = for some c = cq. To calculate the changes in 
the trajectory 5v{u),5s{u) as c is changed by 5c, we write the equations in 
variations for 5v and 5s obtained from Eqs. (|2^ and (p5|) 



d r aiu) ^ 

—5s = 29 

du Vq 

d ^ sl{l~u)^ 3sl{l-u)^ 8s|](1-m)^ 

—5v = „ ' 5v ^ -5s + -5c. (30) 

du C^Vq C^Vo C^Vq 

Since the change in c does not affect the initial conditions, we should have 

5v{0) = 0, 5s{0) = 0. (31) 



Note that according to Eqs. (p8D we have f o ~ m and Sq ~ u, so 5s ~ m and 
5f ~ in the neighborhood of m = 0. 

According to Eq. (0), when u — > we have > for 5c > 0, so 

5v > 0. In turn, according to Eq. (^), -^5s < and therefore 5s < 0. This 
means that the derivatives -^5v and -^5s do not change signs, so 5v > 
everywhere for 5c > and vice versa. Therefore, when the value of c changes 
from to oo, the point at which the trajectory crosses either the line u = 1 
or the line f = in the u — v plane will go monotonically from m = 0, f = 
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Figure 4: The speed c of the travehng pulse as a function of T. The solid 
line shows the results of the numerical solution of the full HH model. The 
dashed line is the prediction of Eq. (p2D. The dotted line is the result of the 
solution of the HH model without the h and n dynamics. 



to M = l,f = 1. Since it depends continuously on c, there is a unique value 
of c = c* at which this point coincides with m = l,w = 0. Numerically, the 
value of c* = I up to the fourth digit. Thus, the dynamics in the pulse front 
uniquely determines its propagation speed. 

Thus, we have obtained an approximate analytical expression for the 
speed of the traveling pulses in the HH model: 

2 / o^o^^nA \ ^'^ 

Equation (|32| ) predicts the speed of the traveling pulse as a function of the 
parameters. To compare this predicted speed with the results obtained from 
the numerical solution of the HH model, we plot the speed c as a function of 
temperature obtained from Eq. (|32[) and from the numerical simulations of 
the HH model (see also ([Huxley, 1959| )) in Fig. § (recall that the temperature 



dependence is contained in the value of am)- As can be seen from this figure, 
the approximate expression for the speed of the pulse given by Eq. (|32D agrees 
with the results for the HH model within 20% at temperatures below 15° C. 
We would like to emphasize that this is the kind of accuracy with which the 
HH model itself predicts the speeds of the traveling pulses as compared to 
the experiments. At higher temperatures the agreement between Eq. (p2D 
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Figure 5: The speed c of the travehng pulse as a function of the membrane 
capacitance C. The sohd hne is the result of the numerical solution of the 
HH model, the dashed line is the prediction of Eq. (|32D. 



and the results of the simulations of the HH model becomes worse, and at 
the temperatures of the propagation threshold Eq. (^) fails completely. We 
have also checked that Eq. predicts the correct dependences on the other 
parameters with similar accuracy at low enough temperatures. For example. 
Fig. ^ shows a comparison of the prediction of Eq. ( |52D with the results 
of the numerical simulations of the HH model at T = 6.3° C as the value 
of the membrane capacitance C per unit area is varied, on the log- log plot. 
Note that the slopes of the two graphs in Fig. ^ agree very well with each 
other. The agreement of the slopes is not as good for the log-log plot of 
the dependence of c on g^^ with other parameters fixed. This is probably 
the consequence of the fact that the errors introdiced by our approximation 

1 /8 

depend on g^g_ stronger than the prediction of the approximation c ~ Qj^^ . 

Incidentally, if the factor of 2/3 in Eq. (^) is replaced by 5/9, it will give 
the the speed of the pulse within a few per cent of that found in the full HH 
model for T < 15° C. Note that if one assumes that m is the fastest variable 
(IFitzHugh, 1961| ; [Casten et al., 19751 ; [Carpenter, 1977| ; [Carpenter, 1979|) and 
calculates the speed of the traveling wave, one will get the value which is an 
order of magnitude greater than the actual value. 

Observe that Fig. ^ also shows the dependence of the speed of the pulse 
on temperature obtained from the simulations of the HH model without the 
h and n dynamics (the dotted line). Note that the insignificance of these 
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variables is one of the key assumptions in deriving Eq. P^). One can see 
that this solution gives an even better approximation to the speed of the 
pulse. This problem, however, cannot be treated analytically in the same 
manner as that for Eqs. ([T5|) and (|l^). 



Let us emphasize that the existence of the front solution is essentially 
determined by the complicated interplay of the V and m kinetics, so the 
problem does not reduce to simple phase plane analysis, in contrast to most 
studies of the traveling waves in excitable systems ([FitzHugh, 1961| ; [Rinzel 



and Keller, 197^ ; [Casten et al., 1975| ; Parpenter, 1977| ; [Carpenter, 1979 



[Vasiliev et al., 1987| ; [Mikhailov, 19901 ). Note that similar situation takes 



place in a class of excitable systems in which the so-called spike autosolitons 
are realized ( |Usipov and IVluratov, 19^5| |Muratov and Osipov, 199UD . These 



models also give rise to the complicated front structures that are similar to 
the one realized in the HH model. 

The validity of the approximations made by us is violated in two cases. 
First, when the temperature becomes sufficiently high, the dynamics of the 
m variable becomes faster, so the separation of the time scales and Ty 
used in our arguments will no longer be justified. One of the implications of 
the absence of this scale separation is the fact that the characteristic value 
of m = m in the front can no longer be treated as small. This allows us to 
derive a criterion for the validity of our approximations 

^ < 1, (33) 

which is equivalent to m < 1 [see Eq. (|2T|)]. In the case of the squid giant 
axon this criterion shows the applicability of our results up to T ^ 25° C, in 
good agreement with Fig. |. 

Another problem may occur when the temperature becomes too low and 
the variable m too slow. In this case the effective time scale r of the dynamics 
in the front of the pulse slows down [see Eq. (pO])], so at some point it may 
become comparable to the leakage time scale r; ~ C/(7; ~ 3 msec. In this 
case one can no longer discard the leakage and the K"*" contributions to the 
membrane current in Eq. ([10|) . Thus, the second validity criterion becomes 
[see Eq. (|g)] 

+ 9K<r < ^ (34) 
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For the squid giant axon, this quantity becomes comparable to 1 only for the 
unrealistically low temperatures T < —30° C. 

As can be seen from Eqs. ( ^31) and (|34D , the quality of the approximations 
used by us should increase with the increase of (^Na- In fact, our procedure 
for finding the spike's speed and the front profile is the leading order of the 
asymptotic limit mg and (^Na ~* oo. 

Figure P shows the functions u(^) and s(^) for c = c* obtained numerically. 
One can see that u{^) has the form of a front connecting the rest state u = 
with the excited state u = 1 which in the original variables corresponds to 
the saturation value V = VNa- 

As can be seen from Fig. ^, the distribution s(^) behind the front ap- 
proaches a straight line with the slope — c*. In terms of the original variables, 
this should correspond to the unlimited growth of m behind the front. This 
however, should not be a problem since this happens only when m ~ m <^ 1. 
When m becomes of order 1, the approximations used to derive Eq. ( PB]) 
ceases to be valid. On the other hand, when this happens, V should already 
be very close to T^ra, so on the time scale Tm <^ T/i.n the variable m will simply 
exponentially relax to m = moo(VNa) behind the front, where, as usual 

This will happen at distances of order cTm S> /, since t [see Eq. (pO|)]. 

If we assume that on the time scale the front was located at z = 0, the 
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Figure 7: The numerical solution of Eqs. (11 



), and (|38|) . 



distribution of m that properly matches with that in Fig. |^ will be 

m{z) = moo(VNa){l - exp[z/crm(VNa)]}, 



(36) 



As was already noted, in the back of the spike and in the refractory tail the 
voltage V changes substantially slower than in the front. Since this happens 
when m ~ 1, the voltage is indeed the fastest variable, so we can eliminate 
it adiabatically from the equations. If we put all the derivatives in Eq. (|10D 
to zero, we find 

gtia.m^hVKi, + QKn^VK + giVi 



V 



(37) 



g^-jn'^h + gvji^ + gi 

This expression uniquely determines the value of ^ as a function of m, h, 
and n. 

To find the approximate distributions of all the variables in the back and 
behind the spike we simply need to solve the initial value problem given by 
Eqs. (|TTD - (|13D with c given by Eq. (|32|) and the following initial conditions 



m(0) = 1, h{Q) = hQ, n(0)=no. 



(38) 



where we assumed that on the even larger length scale CTh^n the front is 
located at 2; = 0. This initial value problem can be straightforwardly solved 
numerically. The result of this solution is shown in Fig. |^ Note that the 
changes in temperature will only change the characteristic length of this 
solution and not its shape. 
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Figure 8: The numerical solution of Eqs. (0), (0), (0), and ( pS] ) with 
m = mooiy) given by Eq. (|35D . 



One can simplify the procedure of finding the distributions of m, h, and 
n by using the fact that -C „ by adiabatically eliminating m. This will 
amount to replacing m by rrioo from Eq. (^) in Eq. (^) and then solving for 
\^ as a function of h and n. The result of the numerical solution of Eqs. (|T^ 
and (p!3D with these approximations is shown in Fig. p. This figure shows 
a good agreement of the slow variables h and n in the refractory tail. Also, 
observe an abrupt jump in the back of the spike. This is due to the fact that 
now the value of V is not uniquely determined by h and n, so at some point 
in the solution a jump occurs from one branch of the dependence V{h, n) to 
the other (see also ([Carpenter, 1977| ; Parpenter, 1979| )). Note that while the 
adiabatic elimination of m works well for the refractory tail, it fails in the 
back of the pulse (compare Figs. ^ to H). 

The results in Figs. |^, |^, and Eq. ( |5B| ) can be combined to give a quanti- 
tative approximation for the whole pulse. This is done in Fig. ^for T = 6.3° 
C. One can see remarkable similarity between the solution of the full HH 
model shown in Fig. || and the one shown in Fig. |^. Thus, our approxima- 
tion scheme has been able to capture quantitatively the essential features of 
the traveling pulses in the HH model. 



4 Conclusions 



We have developed a method that allows to approximately compute the shape 
and the parameters of the traveling spikes in the HH model of conductance 
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Figure 9: An approximate solution for the entire pulse. 



along an axon. Our method is different from the conventional approach 
(IFitzHugh, 19611 ; [Rinzel and Keller, 1973| ; [Casten et al., 1975| ; [Carpenter] 
1977| ; [Carpenter, 1979| ) in the fact that it treats the membrane potential 



rather than the sodium activation variable as the fast variable. We show that 
this is in fact the case for the typical set of the parameters of the Hodgkin- 
Huxley model. This leads to a good quantitative agreement between the 
predictions of the theory and the results of the numerical simulations of the 
HH model. 

Let us emphasize that the HH model itself gives only an approximate, al- 
though quite accurate description of the dynamics of the action potential in 
an actual axon. What we find is that in a broad range of the parameters the 
approximation introduced by us gives an error which is in fact comparable 
to the error produced by the HH model itself as opposed to the experiments. 
For example, at T = 18.5°C the speed of the pulse in the squid giant axon 
was found to be 21.2 m/s ( [Hodgkin and Huxley, 1952 ). The direct numerical 
simulation of the HH model produces the speed of 18.8 m/s, while our proce- 
dure which for this temperature is already near the limit of its applicability 
gives 24.7 m/s. Therefore, we suggest that the ideas of our analysis can in 
fact be built into simpler and more tractable models of nerve conductance 
which will yet be able to give quantitative agreement with the experimental 
observations. 

One of the observations from the analysis made by us is the fact that 
the speed of the traveling spikes in the HH model depends very weakly on 
the slow state variables of the membrane. Indeed, according to Eq. (p2[), the 
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speed of the spike is independent of the vakie of n in front of the spike and 

1 /8 

is proportional to h^' , so a change by a factor of 2 in ho will result in only 
a 10% change in the speed. This makes a perfect biological sense. Thus, 
propagation of the nerve pulses is indeed a very robust phenomenon. 

Another observation one can make from Eq. (|3^) is that if one assumes 
that in addition to the membrane conductance C there is an extra capac- 
itance associated with each sodium channel, there exists a density of the 
channels at which the speed is maximal ( [Hodgkin, 1975|) . Indeed, let us as- 
sume that C = Co + NC^^ and ^^Na = Ng^^^, where Co is the capacitance 
of the membrane without the channels, is the channel density, C^^ is the 
capacitance associated with a single channel, and g^^ is the maximum con- 
ductance of a single channel. For the squid giant axon these parameters are 
estimated to be Co = 0.8 /iF/cm^ C* ^ = 4 x 10"^^ F, c/J^^ = 2.4 x lO'^^ ^-i^ 
and N = 500 /im~^ ( [Hodgkin, 1975| ). Substituting these expressions into 
Eq. (|32D, one gets the speed of the pulse as a function of N. It is not difficult 
to see that this function has a maximum at iV = N^ax = Co/ (4(7^^)- For the 
numerical values above we find N^ax = 500 /im~^, what suggests that the 
channel density in the axon is indeed at the optimum level. The fact that we 
get exactly the same value of N as the one observed may be a bit fortuitous 
because of the approximate nature of Eq. (|3^). Note that because of the very 
slow dependence of the speed on (yfNa, the maximum of the dependence c{N) 
is in fact very fiat, so a change of by a factor of 2 from N^ax results only 
in a 7% difference in c. 

So far, we were talking only about the traveling wave solutions in the 
form of the solitary spikes. It is not difficult to see that our method can be 
extended to spike trains as well. Indeed, the speed of a spike in a spike train 
is determined by the value of the slow variable h in front of the spike [see 
Eq. (P^j, which, however, is now different from the equilibrium value ho and 
must be determined. Outside of the spike fronts one has to solve the equations 
of the slow dynamics given by Eqs. ( |Tl| ) - (|T3|) in which the fast variable V 
has been eliminated adiabatically via Eq. (pT]). These equations have to be 
solved as an initial value problem for —L < z < with m(0) = rriooiVj^sg), 
h{0) = hg, and n(0) = Ug. Here hg and rig are the values of h and n in the 
spike, respectively, L is the spatial period of the spike train, and we assumed 
that the front of one of the spikes in the spike train is located at z = 0. The 
spikes are also assumed to travel to the right with the speed given by Eq. (|32D 
in which ho is replaced by hg. Then, the values of hg and have to be found 
self-consistently from the condition that h{—L) = hg and n{—L) = Hg. 



19 



We implemented this procedure numerically. To find the values of hg and 
Hs as functions of L, we started with a reasonable initial guess for hg and 
and then solved the initial value problem up to z = —L. The values of h and 
n at z = —L were then taken as the new values for hg and Ug, respectively, 
and the procedure was iterated. We found a fast convergence to the periodic 
solution. Then, from the value of hg we obtained the speed of the spike 
train as a function of L and therefore the dependence of the spike frequency 
on the period. Our main finding was that in the region of the applicability 
of our approximation [Eq. (|3^) ] the speed of the spike trains is practically 
independent of the period, so they have almost no dispersion. This also 
makes good biological sense. In fact, the amount of the dispersion we found 
is comparable to the error introduced by our approximation scheme. Since 
we are only interested in quantitative predictions, we do not present these 
results in detail here. Also, our method fails for periods L < 10 cm, for which 
a substantial amount of dispersion was found in the simulations of the full 
HH model ( [Miller and Rinzel, 1981|) . Nevertheless, let us point out that the 
results obtained with our method are in a good qualitative agreement with 
those obtained for the full HH model ( [Miller and Rinzel, 1981|) . In particular, 
according to our numerical solution outlined above there exists a period of 
the spike train for which the speed of the spikes reaches maximum, greater 
than that of the solitary spike due to a slight overshoot of the h variable 
behind the spike. However, the magnitude of this overshoot is so small that 
it only changes the speed of the spike by a fraction of a per cent. So, for 
practical purposes the spike trains with period L > 10 cm may be considered 
dispersionless. 

In short, we have introduced an approximation scheme which allows to 
make quantitative predictions of the shape and the parameters of the travel- 
ing pulses in the HH model of nerve conductance. We hope that our results 
will provide an easy and convenient tool for analyzing the fascinating com- 
plexity of the neural activity. 
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